      SUBROUTINE  DEIGQR ( A,N,NA,H,NH,IPARAM,MULT,INT,LAMBDA,IER )     00010001
C                                                                       00011001
CCCCC --------EIGENVALUES ANG EIGENFUNCTIONS FOR NON-HERMITIAN,---------00020000
CCCCC         COMPLEX MATRIXES BY THE QR-METHOD.                        00030000
CCCCC ---------------- COPIED FROM SYS1.MSL(DEIGCN)  ------------------ 00040000
C                                                                       00041001
      IMPLICIT REAL*8 (A-H,O-Z)                                         00050000
      COMPLEX*16 A,LAMBDA,MULT,H,SFT(2),SHIFT,TEMP,TEMP1,TEMP2,         00060000
     2           ASIN,ACOS,BSIN,BCOS,ASUM,AMUL,ZN                       00070000
      LOGICAL TWICE                                                     00080000
      DIMENSION A(NA,1),LAMBDA(1),MULT(1),H(NH,1),INT(1)                00090000
C                                                                       00100000
C     INITIALIZE                                                        00110000
C                                                                       00120000
CC    CALL ERRSET(207,256,-1,1)                                         00130000
CC     CALL ERRSET(208,256,-1,1)                                        00140000
      IF( N.LE.0 ) GO TO 70                                             00150000
      IF( NA.LT.N ) GO TO 70                                            00160000
       IF( NH.LT.NA*2 ) GO TO 70                                        00170000
      IER=0                                                             00180000
      M=N                                                               00190000
      IF(N.NE.1) GO TO 1                                                00200000
      LAMBDA(1)=A(1,1)                                                  00210000
      A(1,1)=(1.D0,0.D0)                                                00220000
      RETURN                                                            00230000
C                                                                       00240000
    1 ICOUNT=0                                                          00250000
      N1=N-1                                                            00260000
      NCAL=N                                                            00270000
      SHIFT=(0.D0,0.D0)                                                 00280000
      IF(N.NE.2) GO TO 4                                                00290000
C                                                                       00300000
C     EIGENVALUES OF 2*2 MATRIX                                         00310000
C                                                                       00320000
    2 ASUM=A(1,1)+A(2,2)                                                00330000
      AMUL=A(1,1)*A(2,2)-A(1,2)*A(2,1)                                  00340000
      TEMP=(ASUM+CDSQRT(ASUM*ASUM-4.D0*AMUL))*0.5D0                     00350000
      IF( DREAL(TEMP).NE.0.D0 .OR. DIMAG(TEMP).NE.0.D0 ) GO TO 3        00360000
      LAMBDA(M)=SHIFT                                                   00370000
      LAMBDA(M-1)=ASUM+SHIFT                                            00380000
      GO TO 37                                                          00390000
C                                                                       00400000
C     REDUCE ORIGINAL MATRIX TO HESSENBERG FORM                         00410000
C                                                                       00420000
    4 N2=N-2                                                            00430000
      DO 15 K=1,N2                                                      00440000
      K1=K+1                                                            00450000
      K2=K+2                                                            00460000
C     FIND THE LARGEST ELEMENT IN EACH COLUMN                           00470000
      ABIG=0.D0                                                         00480000
      INT(K)=K1                                                         00490000
      DO 5 I=K1,N                                                       00500000
      ABSSQ=DREAL(A(I,K))**2+DIMAG(A(I,K))**2                           00510000
      IF(ABSSQ.LE.ABIG) GO TO 5                                         00520000
      INT(K)=I                                                          00530000
      ABIG=ABSSQ                                                        00540000
    5 CONTINUE                                                          00550000
C     EXCHANGE ROWS IF NECESSARY                                        00560000
      INTER=INT(K)                                                      00570000
      IF(ABIG.EQ.0.D0) GO TO 15                                         00580000
      IF(INTER.EQ.K1) GO TO 8                                           00590000
      DO 6 J=K,N                                                        00600000
      TEMP=A(K1,J)                                                      00610000
      A(K1,J)=A(INTER,J)                                                00620000
      A(INTER,J)=TEMP                                                   00630000
    6 CONTINUE                                                          00640000
C     EXCHANGE CORRESPONDING COLUMNS                                    00650000
      DO 7 I=1,N                                                        00660000
      TEMP=A(I,K1)                                                      00670000
      A(I,K1)=A(I,INTER)                                                00680000
      A(I,INTER)=TEMP                                                   00690000
    7 CONTINUE                                                          00700000
C     ELIMINATE                                                         00710000
    8 ZN=1.D0/A(K1,K)                                                   00720000
      DO 9 I=K2,N                                                       00730000
      MULT(I)=A(I,K)*ZN                                                 00740000
      A(I,K)=MULT(I)                                                    00750000
    9 CONTINUE                                                          00760000
      DO 11 J=K2,N                                                      00770000
      DO 10 I=1,K1                                                      00780000
      A(I,K1)=A(I,K1)+A(I,J)*MULT(J)                                    00790000
   10 CONTINUE                                                          00800000
   11 CONTINUE                                                          00810000
      DO 13 J=K2,N                                                      00820000
      DO 12 I=K2,N                                                      00830000
      A(I,K1)=A(I,K1)+A(I,J)*MULT(J)                                    00840000
   12 CONTINUE                                                          00850000
   13 CONTINUE                                                          00860000
      DO 60 I=K2,N                                                      00870000
      A(I,K1)=A(I,K1)-MULT(I)*A(K1,K1)                                  00880000
   60 CONTINUE                                                          00890000
      DO 14 J=K2,N                                                      00900000
      DO 14 I=K2,N                                                      00910000
      A(I,J)=A(I,J)-MULT(I)*A(K1,J)                                     00920000
   14 CONTINUE                                                          00930000
   15 CONTINUE                                                          00940000
C                                                                       00950000
C     CALCULATE EPSILON                                                 00960000
C                                                                       00970000
      EPS=0.D0                                                          00980000
      DO 16 J=1,N                                                       00990000
      EPS=EPS+CDABS(A(1,J))                                             01000000
   16 CONTINUE                                                          01010000
      DO 18 I=2,N                                                       01020000
      SUM=0.D0                                                          01030000
      I1=I-1                                                            01040000
      DO 17 J=I1,N                                                      01050000
      SUM=SUM+CDABS(A(I,J))                                             01060000
   17 CONTINUE                                                          01070000
      IF(SUM.GT.EPS) EPS=SUM                                            01080000
   18 CONTINUE                                                          01090000
      E=DSQRT(DFLOAT(N))                                                01100000
      EPS=E*EPS*1.0D-20                                                 01110000
      IF(EPS.EQ.0.D0) EPS=1.D-20                                        01120000
      IF(IPARAM.EQ.0) GO TO 20                                          01130000
C     SAVE ELEMENTS OF HESSENBERG MATRIX                                01140000
      DO 19 J=1,N                                                       01150000
      DO 19 I=1,N                                                       01160000
      H(N+I,J)=A(I,J)                                                   01170000
   19 CONTINUE                                                          01180000
C                                                                       01190000
C     START ITERATION                                                   01200000
C                                                                       01210000
   20 IF(N.NE.1) GO TO 21                                               01220000
      LAMBDA(M)=A(1,1)+SHIFT                                            01230000
      GO TO 37                                                          01240000
   21 IF(N.EQ.2) GO TO 2                                                01250000
   22 IF(DREAL(A(N,N)).EQ.0.D0 .AND. DIMAG(A(N,N)).EQ.0.D0 ) GO TO 23   01260000
      TEMP=A(N,N1)/A(N,N)                                               01270000
      IF(DABS(DREAL(TEMP))+DABS(DIMAG(TEMP))-1.D-16) 24,24,23           01280000
   23 IF(DABS(DREAL(A(N,N1)))+DABS(DIMAG(A(N,N1))).GE.EPS) GO TO 25     01290000
   24 LAMBDA(M-N+1)=A(N,N)+SHIFT                                        01300000
      ICOUNT=0                                                          01310000
      N=N-1                                                             01320000
      N1=N-1                                                            01330000
      GO TO 21                                                          01340000
C                                                                       01350000
C     DETERMINE SHIFT                                                   01360000
C                                                                       01370000
   25 ASUM=A(N1,N1)+A(N,N)                                              01380000
      AMUL=A(N1,N1)*A(N,N)-A(N1,N)*A(N,N1)                              01390000
      TEMP=ASUM*ASUM-4.D0*AMUL                                          01400000
      SFT(1)=(ASUM+CDSQRT(TEMP))*0.5D0                                  01410000
      SFT(2)=AMUL/SFT(1)                                                01420000
      IF(DREAL(SFT(1)).EQ.0.D0 .AND. DIMAG(SFT(1)).EQ.0.D0 ) SFT(2)=ASUM01430000
      TEMP1=SFT(1)-A(N,N)                                               01440000
      TEMP2=SFT(2)-A(N,N)                                               01450000
      INDEX=1                                                           01460000
      IF(CDABS(TEMP1).GE.CDABS(TEMP2)) INDEX=2                          01470000
      IF(CDABS(A(N1,N-2)).GE.EPS) GO TO 30                              01480000
      LAMBDA(M-N+1)=SFT(1)+SHIFT                                        01490000
      LAMBDA(M-N+2)=SFT(2)+SHIFT                                        01500000
      ICOUNT=0                                                          01510000
      N=N-2                                                             01520000
      N1=N-1                                                            01530000
      GO TO 20                                                          01540000
   30 SHIFT=SHIFT+SFT(INDEX)                                            01550000
      DO 31 I=1,N                                                       01560000
      A(I,I)=A(I,I)-SFT(INDEX)                                          01570000
   31 CONTINUE                                                          01580000
C                                                                       01590000
C     PERFORM GIVENS ROTATIONS, QR ITERATIONS                           01600000
C                                                                       01610000
      IF(ICOUNT.LE.10) GO TO 32                                         01620000
      NCAL=M-N                                                          01630000
      IER=1                                                             01640000
      GO TO 37                                                          01650000
   32 DO 36 K=1,N1                                                      01660000
      K1=K+1                                                            01670000
      K2=K+2                                                            01680000
      TEMP1=A(K,K)                                                      01690000
      TEMP2=A(K1,K)                                                     01700000
      ROOT=DSQRT(DREAL(TEMP1)**2+DIMAG(TEMP1)**2                        01710000
     1          +DREAL(TEMP2)**2+DIMAG(TEMP2)**2)                       01720000
      IF(ROOT.EQ.0.D0) GO TO 36                                         01730000
      ACOS=TEMP1/ROOT                                                   01740000
      ASIN=TEMP2/ROOT                                                   01750000
      BCOS=DCONJG(ACOS)                                                 01760000
      BSIN=DCONJG(ASIN)                                                 01770000
      INDEX=MAX0(K-1,1)                                                 01780000
      DO 33 J=INDEX,N                                                   01790000
      TEMP=A(K,J)                                                       01800000
      A(K,J)=BCOS*A(K,J)+BSIN*A(K1,J)                                   01810000
      A(K1,J)=-ASIN*TEMP+ACOS*A(K1,J)                                   01820000
   33 CONTINUE                                                          01830000
      DO 34 I=1,K                                                       01840000
      TEMP=A(I,K)                                                       01850000
      A(I,K)=ACOS*A(I,K)+ASIN*A(I,K1)                                   01860000
      A(I,K1)=-BSIN*TEMP+BCOS*A(I,K1)                                   01870000
   34 CONTINUE                                                          01880000
      INDEX=MIN0(K2,N)                                                  01890000
      DO 35 I=K1,INDEX                                                  01900000
      A(I,K)=ASIN*A(I,K1)                                               01910000
      A(I,K1)=BCOS*A(I,K1)                                              01920000
   35 CONTINUE                                                          01930000
   36 CONTINUE                                                          01940000
      ICOUNT=ICOUNT+1                                                   01950000
      GO TO 22                                                          01960000
    3 LAMBDA(M)=TEMP+SHIFT                                              01970000
      LAMBDA(M-1)=AMUL/TEMP+SHIFT                                       01980000
   37 N=M                                                               01990000
      IF(NCAL.EQ.0) GO TO 59                                            02000000
      IF(IPARAM.EQ.0) RETURN                                            02010000
C                                                                       02020000
C     CALCULATE VECTORS                                                 02030000
C                                                                       02040000
      N1=N-1                                                            02050000
      IF(N.NE.2) GO TO 38                                               02060000
      EPS=DMAX1(CDABS(LAMBDA(1)),CDABS(LAMBDA(2)))*1.D-16               02070000
      IF(EPS.EQ.0.D0) EPS=1.D-20                                        02080000
      H(3,1)=A(1,1)                                                     02090000
      H(3,2)=A(1,2)                                                     02100000
      H(4,1)=A(2,1)                                                     02110000
      H(4,2)=A(2,2)                                                     02120000
   38 DO 57 K=1,NCAL                                                    02130000
      DO 40 J=1,N                                                       02140000
      DO 39 I=1,N                                                       02150000
      H(I,J)=H(N+I,J)                                                   02160000
   39 CONTINUE                                                          02170000
   40 CONTINUE                                                          02180000
      DO 61 I=1,N                                                       02190000
      H(I,I)=H(I,I)-LAMBDA(K)                                           02200000
   61 CONTINUE                                                          02210000
C     REDUCE HESSENBERG MATRIX TO UPPER TRIANGULAR FORM                 02220000
      DO 44 I=1,N1                                                      02230000
      I1=I+1                                                            02240000
      MULT(I)=(0.D0,0.D0)                                               02250000
      INT(N+I)=0                                                        02260000
      IF(CDABS(H(I1,I)).LE.CDABS(H(I,I))) GO TO 42                      02270000
      INT(N+I)=1                                                        02280000
      DO 41 J=I,N                                                       02290000
      TEMP=H(I,J)                                                       02300000
      H(I,J)=H(I1,J)                                                    02310000
      H(I1,J)=TEMP                                                      02320000
   41 CONTINUE                                                          02330000
   42 IF(DREAL(H(I,I)).EQ.0.D0 .AND. DIMAG(H(I,I)).EQ.0.D0) GO TO 44    02340000
      MULT(I)=-H(I1,I)/H(I,I)                                           02350000
      DO 43 J=I1,N                                                      02360000
      H(I1,J)=H(I1,J)+MULT(I)*H(I,J)                                    02370000
   43 CONTINUE                                                          02380000
   44 CONTINUE                                                          02390000
C     START INVERSE ITERATION                                           02400000
      DO 45 I=1,N                                                       02410000
      IF(DREAL(H(I,I)).EQ.0.D0 .AND. DIMAG(H(I,I)).EQ.0.D0)             02420000
     1   H(I,I)=DCMPLX(EPS,0.D0)                                        02430000
      A(I,K)=(1.D0,0.D0)                                                02440000
   45 CONTINUE                                                          02450000
      TWICE=.FALSE.                                                     02460000
   46 A(N,K)=A(N,K)/H(N,N)                                              02470000
      DO 48 I=1,N1                                                      02480000
      NI=N-I                                                            02490000
      DO 47 J=NI,N1                                                     02500000
      A(NI,K)=A(NI,K)-H(NI,J+1)*A(J+1,K)                                02510000
   47 CONTINUE                                                          02520000
      A(NI,K)=A(NI,K)/H(NI,NI)                                          02530000
   48 CONTINUE                                                          02540000
      IF(TWICE) GO TO 52                                                02550000
      BIG=0.D0                                                          02560000
      DO 49 I=1,N                                                       02570000
      SUM=DABS(DREAL(A(I,K)))+DABS(DIMAG(A(I,K)))                       02580000
      IF(BIG.LT.SUM) BIG=SUM                                            02590000
   49 CONTINUE                                                          02600000
      ZN=1.D0/BIG                                                       02610000
      DO 50 I=1,N                                                       02620000
      A(I,K)=A(I,K)*ZN                                                  02630000
   50 CONTINUE                                                          02640000
      DO 51 I=1,N1                                                      02650000
      IF(INT(N+I).EQ.0) GO TO 51                                        02660000
      TEMP=A(I,K)                                                       02670000
      A(I,K)   = A(I+1,K)                                               02680000
      A(I+1,K) = TEMP                                                   02690000
   51 A(I+1,K) = A(I+1,K)+MULT(I)*A(I,K)                                02700000
      TWICE=.TRUE.                                                      02710000
      GO TO 46                                                          02720000
C     TRANSFORM EIGENVECTORS                                            02730000
   52 IF(N.EQ.2) GO TO 55                                               02740000
      N2=N-2                                                            02750000
      DO 54 J=1,N2                                                      02760000
      NP=N-J+1                                                          02770000
      DO 53 I=NP,N                                                      02780000
      A(I,K) = A(I,K)+H(N+I,N-J-1)*A(N-J,K)                             02790000
   53 CONTINUE                                                          02800000
      INDEX    = INT(N-J-1)                                             02810000
      TEMP     = A(N-J,K)                                               02820000
      A(N-J,K) = A(INDEX,K)                                             02830000
      A(INDEX,K)=TEMP                                                   02840000
   54 CONTINUE                                                          02850000
C     NORMALIZE                                                         02860000
   55 BIG=0.D0                                                          02870000
      DO 56 I=1,N                                                       02880000
      SUM=CDABS(A(I,K))                                                 02890000
      IF(BIG.GE.SUM) GO TO 56                                           02900000
      BIG=SUM                                                           02910000
      TEMP=A(I,K)                                                       02920000
   56 CONTINUE                                                          02930000
      ZN=1.D0/TEMP                                                      02940000
      DO 57 I=1,N                                                       02950000
      A(I,K)=A(I,K)*ZN                                                  02960000
   57 CONTINUE                                                          02970000
      IF(IER.EQ.0) RETURN                                               02980000
      NCAL1=NCAL+1                                                      02990000
      DO 58 K=NCAL1,N                                                   03000000
      DO 58 I=1,N                                                       03010000
      A(I,K)=(0.D0,0.D0)                                                03020000
   58 CONTINUE                                                          03030000
      RETURN                                                            03040000
   59 IER=2                                                             03050000
      RETURN                                                            03060000
   70 IER=3                                                             03070000
      WRITE(6,1000) N,NA,NH                                             03080000
 1000 FORMAT(1H0,'(SUBR. DEIGCN) INVALID ARGUMENT. N,NA,NH = ',I5,      03090000
     2           ',',I5,',',I5)                                         03100000
      RETURN                                                            03110000
      END                                                               03120000
